Lecture 14: Spatial GLMs

Class Intro

Intro Questions

  • Describe geometric anisotropy. Discuss this in the context of a data set.

  • For Today:
    • Spatial Generalized Linear Models

Geometric Anisotropy

Geometric Anisotropy Model

  • Let \(\boldsymbol{Y}(\boldsymbol{s}) = \mu(\boldsymbol{s}) + w(\boldsymbol{s}) + \epsilon(\boldsymbol{s})\), thus \(\boldsymbol{Y}(\boldsymbol{s}) \sim N(\mu(\boldsymbol{s}), \Sigma(\tau^2, \sigma^2, \phi, B))\), where \(B = L^T L\).
  • The covariance matrix is defined as \(\Sigma(\tau^2, \sigma^2, \phi, B)) = \tau^2 I + \sigma^2 H((\boldsymbol{h}^T B \boldsymbol{h}^T)^{\frac{1}{2}}),\) where \(H((\boldsymbol{h}^T B \boldsymbol{h}^T)^{\frac{1}{2}})\) has entries of \(\rho((\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}}))\) with \(\rho()\) being a valid correation function, typically including \(\phi\) and \(\boldsymbol{h_{ij}} = \boldsymbol{s_i} - \boldsymbol{s_j}\).

Geometric Anisotropy Visual

  • Consider four points positioned on a unit circle.

  • How far apart are each set of points?

Geometric Anisotropy Exercise 1

Now consider a set of correlation functions. For each, calculate the correlation matrix and discuss the impact of \(B\) on the correlation. Furthermore, how does B change the geometry of the correlation?

  1. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

  2. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 2 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

  3. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 3 & 1 \\ 1 & 1 \\ \end{pmatrix}\)

Geometric Anisotropy: Solution 1

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T I \boldsymbol{h_{ij}}^T)^{\frac{1}{2}}))\)

Implied Distance

0.00 1.41 1.41 2.00
1.41 0.00 2.00 1.41
1.41 2.00 0.00 1.41
2.00 1.41 1.41 0.00

Correlation

1.000 0.243 0.243 0.135
0.243 1.000 0.135 0.243
0.243 0.135 1.000 0.243
0.135 0.243 0.243 1.000

Geometric Anisotropy: Solution 2

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 2 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

Implied Distance

0.00 1.73 1.73 2.83
1.73 0.00 2.00 1.73
1.73 2.00 0.00 1.73
2.83 1.73 1.73 0.00

Correlation

1.000 0.177 0.177 0.059
0.177 1.000 0.135 0.177
0.177 0.135 1.000 0.177
0.059 0.177 0.177 1.000

Geometric Anisotrop: Solution 3

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 3 & 1 \\ 1 & 1 \\ \end{pmatrix}\)

Implied Distance

0.00 1.41 2.45 3.46
1.41 0.00 2.00 2.45
2.45 2.00 0.00 1.41
3.46 2.45 1.41 0.00

Correlation

1.000 0.243 0.086 0.031
0.243 1.000 0.135 0.086
0.086 0.135 1.000 0.243
0.031 0.086 0.243 1.000

More Geometric Anisotropy

  • The matrix \(B\) relates to the orientation of a transformed ellipse.
  • The (effective) range for any angle \(\eta\) is determined by the equation \[\rho(r_\eta(\tilde{\boldsymbol{h}}_{\eta}^T B \tilde{\boldsymbol{h}}_{\eta}^T)^{\frac{1}{2}}) = .05,\] where \(\tilde{\boldsymbol{h}}_{\eta}\) is a unit vector in the direction \(\eta\).

Fitting Geometric Anisotropy Models

  • Okay, so if we suspect that geometric anisotrophy is present, how do we fit the model? That is, what is necessary in estimating this model?
  • In addition to \(\sigma^2\) and \(\tau^2\) we need to fit \(B\).
  • What about \(\phi\)? What is \(\phi\) when \(B = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}\)?

Priors for B

  • While \(B\) is a matrix, it is just another unknown parameter.
  • Hence, to fit a Bayesian model we need a prior distribution for \(B\).
  • One option for the positive definite matrix is the Wishart distribution, which is a bit like a matrix-variate gamma distribution.

Spatial GLMS

Data Viz Motivation

Ozone Exceedance in Colorado

Generalized Linear Model Notation

There are three components to a generalized linear model:

  1. Sampling Distribution: such as Poisson or Binomial
  2. Linear combination of predictors: \(\eta = X\beta\)
  3. A link function to map the linear combination of predictors to the support of the sampling distribution.

Logistic Regression Overview

Write out the complete model specification for logistic regression.

  • Assume \(Y_i\) is the binary response for the \(i^{th}\) observation, \[\begin{eqnarray*} Y_i &\sim& Bernoulli(\pi_i)\\ logit(\pi_i) &=& X_i \beta, \end{eqnarray*}\]

  • where \(logit(\pi_i) = log \left(\frac{\pi_i}{1-\pi_i}\right)\)

glm()

Interpret this output

## 
## Call:
## glm(formula = Exceedance ~ north, family = binomial(link = "logit"), 
##     data = CO)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.313   0.378   0.378   0.378   1.177  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -1.695e-15  8.165e-01   0.000   1.0000  
## north        2.603e+00  1.097e+00   2.372   0.0177 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.708  on 34  degrees of freedom
## Residual deviance: 22.873  on 33  degrees of freedom
## AIC: 26.873
## 
## Number of Fisher Scoring iterations: 5

Logistic Regression in JAGS

## 
## Iterations = 2001:7000
## Thinning interval = 1 
## Number of chains = 2 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean      SD  Naive SE Time-series SE
## int   -0.005624 0.88044 0.0088044      0.0207100
## north  2.865766 1.20013 0.0120013      0.0275907
## p[1]   0.931508 0.04541 0.0004541      0.0005299
## p[2]   0.931508 0.04541 0.0004541      0.0005299
## p[3]   0.498883 0.18685 0.0018685      0.0042908
## p[4]   0.931508 0.04541 0.0004541      0.0005299
## p[5]   0.498883 0.18685 0.0018685      0.0042908
## p[6]   0.931508 0.04541 0.0004541      0.0005299
## p[7]   0.931508 0.04541 0.0004541      0.0005299
## p[8]   0.498883 0.18685 0.0018685      0.0042908
## p[9]   0.931508 0.04541 0.0004541      0.0005299
## p[10]  0.931508 0.04541 0.0004541      0.0005299
## p[11]  0.498883 0.18685 0.0018685      0.0042908
## p[12]  0.931508 0.04541 0.0004541      0.0005299
## p[13]  0.498883 0.18685 0.0018685      0.0042908
## p[14]  0.498883 0.18685 0.0018685      0.0042908
## p[15]  0.931508 0.04541 0.0004541      0.0005299
## p[16]  0.931508 0.04541 0.0004541      0.0005299
## p[17]  0.931508 0.04541 0.0004541      0.0005299
## p[18]  0.931508 0.04541 0.0004541      0.0005299
## p[19]  0.931508 0.04541 0.0004541      0.0005299
## p[20]  0.931508 0.04541 0.0004541      0.0005299
## p[21]  0.931508 0.04541 0.0004541      0.0005299
## p[22]  0.931508 0.04541 0.0004541      0.0005299
## p[23]  0.931508 0.04541 0.0004541      0.0005299
## p[24]  0.931508 0.04541 0.0004541      0.0005299
## p[25]  0.931508 0.04541 0.0004541      0.0005299
## p[26]  0.931508 0.04541 0.0004541      0.0005299
## p[27]  0.931508 0.04541 0.0004541      0.0005299
## p[28]  0.931508 0.04541 0.0004541      0.0005299
## p[29]  0.931508 0.04541 0.0004541      0.0005299
## p[30]  0.931508 0.04541 0.0004541      0.0005299
## p[31]  0.931508 0.04541 0.0004541      0.0005299
## p[32]  0.931508 0.04541 0.0004541      0.0005299
## p[33]  0.931508 0.04541 0.0004541      0.0005299
## p[34]  0.931508 0.04541 0.0004541      0.0005299
## p[35]  0.931508 0.04541 0.0004541      0.0005299
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%      50%    75%  97.5%
## int   -1.7251 -0.5796 0.004706 0.5623 1.7255
## north  0.6192  2.0478 2.840331 3.6219 5.3305
## p[1]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[2]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[3]   0.1512  0.3590 0.501176 0.6370 0.8488
## p[4]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[5]   0.1512  0.3590 0.501176 0.6370 0.8488
## p[6]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[7]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[8]   0.1512  0.3590 0.501176 0.6370 0.8488
## p[9]   0.8198  0.9070 0.940798 0.9659 0.9912
## p[10]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[11]  0.1512  0.3590 0.501176 0.6370 0.8488
## p[12]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[13]  0.1512  0.3590 0.501176 0.6370 0.8488
## p[14]  0.1512  0.3590 0.501176 0.6370 0.8488
## p[15]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[16]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[17]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[18]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[19]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[20]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[21]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[22]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[23]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[24]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[25]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[26]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[27]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[28]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[29]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[30]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[31]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[32]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[33]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[34]  0.8198  0.9070 0.940798 0.9659 0.9912
## p[35]  0.8198  0.9070 0.940798 0.9659 0.9912
##                  2.5 %   97.5 %
## (Intercept) -1.6870118 1.687012
## north        0.5038285 4.968462

Spatial Logistic Regression

  • Assume \(Y_i\) is the binary response for the \(i^{th}\) observation, \[\begin{eqnarray*} Y_i|\beta, w(\boldsymbol{s_i}) &\sim& Bernoulli(\pi_i)\\ logit(\pi_i) &=& X_i \beta + w(\boldsymbol{s_i}), \\ \end{eqnarray*}\]
  • where \(\boldsymbol{W} \sim N(\boldsymbol{0},\sigma^2 H(\phi))\)

Spatial Logistic Regression in JAGS

Spatial Logistic Regression with spGLM()

## ----------------------------------------
##  General model description
## ----------------------------------------
## Model fit with 64 observations.
## 
## Number of covariates 1 (including intercept if specified).
## 
## Using the exponential spatial correlation model.
## 
## Number of MCMC samples 10000.
## 
## Priors and hyperpriors:
##  beta normal:
##      mu: 0.000   
##      sd: 10.000  
## 
## 
##  sigma.sq IG hyperpriors shape=2.00000 and scale=1.00000
## 
##  phi Unif hyperpriors a=0.03000 and b=0.30000
## 
## Adaptive Metropolis with target acceptance rate: 43.0
## -------------------------------------------------
##      Sampling
## -------------------------------------------------
## Batch: 10 of 200, 5.00%
##  parameter   acceptance  tuning
##  beta[0]     64.0        0.11900
##  sigma.sq    40.0        0.46620
##  phi     70.0        0.55814
## -------------------------------------------------
## Batch: 20 of 200, 10.00%
##  parameter   acceptance  tuning
##  beta[0]     66.0        0.13151
##  sigma.sq    36.0        0.45697
##  phi     72.0        0.61684
## -------------------------------------------------
## Batch: 30 of 200, 15.00%
##  parameter   acceptance  tuning
##  beta[0]     52.0        0.14534
##  sigma.sq    32.0        0.46620
##  phi     64.0        0.68171
## -------------------------------------------------
## Batch: 40 of 200, 20.00%
##  parameter   acceptance  tuning
##  beta[0]     52.0        0.16063
##  sigma.sq    40.0        0.43905
##  phi     70.0        0.75341
## -------------------------------------------------
## Batch: 50 of 200, 25.00%
##  parameter   acceptance  tuning
##  beta[0]     62.0        0.17752
##  sigma.sq    40.0        0.44792
##  phi     60.0        0.83265
## -------------------------------------------------
## Batch: 60 of 200, 30.00%
##  parameter   acceptance  tuning
##  beta[0]     46.0        0.19619
##  sigma.sq    38.0        0.43035
##  phi     76.0        0.92022
## -------------------------------------------------
## Batch: 70 of 200, 35.00%
##  parameter   acceptance  tuning
##  beta[0]     48.0        0.21683
##  sigma.sq    36.0        0.42183
##  phi     52.0        1.01700
## -------------------------------------------------
## Batch: 80 of 200, 40.00%
##  parameter   acceptance  tuning
##  beta[0]     44.0        0.23489
##  sigma.sq    32.0        0.43905
##  phi     46.0        1.12395
## -------------------------------------------------
## Batch: 90 of 200, 45.00%
##  parameter   acceptance  tuning
##  beta[0]     50.0        0.24941
##  sigma.sq    52.0        0.43905
##  phi     56.0        1.16982
## -------------------------------------------------
## Batch: 100 of 200, 50.00%
##  parameter   acceptance  tuning
##  beta[0]     58.0        0.27018
##  sigma.sq    42.0        0.43035
##  phi     42.0        1.21756
## -------------------------------------------------
## Batch: 110 of 200, 55.00%
##  parameter   acceptance  tuning
##  beta[0]     40.0        0.28121
##  sigma.sq    48.0        0.43905
##  phi     40.0        1.29285
## -------------------------------------------------
## Batch: 120 of 200, 60.00%
##  parameter   acceptance  tuning
##  beta[0]     56.0        0.29860
##  sigma.sq    38.0        0.45697
##  phi     52.0        1.31897
## -------------------------------------------------
## Batch: 130 of 200, 65.00%
##  parameter   acceptance  tuning
##  beta[0]     36.0        0.29860
##  sigma.sq    38.0        0.43905
##  phi     50.0        1.42883
## -------------------------------------------------
## Batch: 140 of 200, 70.00%
##  parameter   acceptance  tuning
##  beta[0]     44.0        0.30463
##  sigma.sq    44.0        0.42183
##  phi     58.0        1.45769
## -------------------------------------------------
## Batch: 150 of 200, 75.00%
##  parameter   acceptance  tuning
##  beta[0]     44.0        0.29860
##  sigma.sq    48.0        0.43905
##  phi     56.0        1.54783
## -------------------------------------------------
## Batch: 160 of 200, 80.00%
##  parameter   acceptance  tuning
##  beta[0]     32.0        0.30463
##  sigma.sq    36.0        0.41348
##  phi     38.0        1.51718
## -------------------------------------------------
## Batch: 170 of 200, 85.00%
##  parameter   acceptance  tuning
##  beta[0]     44.0        0.30463
##  sigma.sq    30.0        0.40529
##  phi     38.0        1.61100
## -------------------------------------------------
## Batch: 180 of 200, 90.00%
##  parameter   acceptance  tuning
##  beta[0]     40.0        0.29860
##  sigma.sq    42.0        0.41348
##  phi     70.0        1.61100
## -------------------------------------------------
## Batch: 190 of 200, 95.00%
##  parameter   acceptance  tuning
##  beta[0]     24.0        0.31079
##  sigma.sq    58.0        0.41348
##  phi     48.0        1.54783
## -------------------------------------------------
## Sampled: 10000 of 10000, 100.00%
## -------------------------------------------------
## 
## Iterations = 9000:10000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 1001 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                Mean     SD Naive SE Time-series SE
## (Intercept) 0.30659 0.2909 0.009195       0.056653
## sigma.sq    1.60258 0.5192 0.016411       0.077912
## phi         0.08361 0.0423 0.001337       0.006326
## 
## 2. Quantiles for each variable:
## 
##                 2.5%     25%     50%     75%  97.5%
## (Intercept) -0.28691 0.12352 0.32349 0.51480 0.8205
## sigma.sq     0.85714 1.21177 1.52577 1.86219 2.8379
## phi          0.03291 0.05516 0.07694 0.09632 0.2161

Poisson Regression Overview

Poisson Regression in JAGS

Spatial Poisson Regression in JAGS

Spatial Poisson Regression with spGLM()